#include <assert.h>
#include <iostream>
#include <numeric>

#include "foreach.hpp"
#include "heightmap.hpp"

namespace hex {

heightmap::heightmap() : width_(0), height_(0), subdivisions_(0), points_between_(0)
{}

heightmap::heightmap(int width, int height)
  : map_(width*height), width_(width), height_(height), subdivisions_(0), points_between_(0)
{}

void heightmap::set(int x, int y, GLfloat value)
{
	assert(x >= 0 && x < width_);
	assert(y >= 0 && y < height_);
	mutable_val(x, y) = value;
}

void heightmap::add(int x, int y, GLfloat value)
{
	assert(x >= 0 && x < width_);
	assert(y >= 0 && y < height_);
	mutable_val(x, y) += value;
}

namespace {

GLfloat calculate_variance(GLfloat variance)
{
	int int_variance = static_cast<int>(variance*1000.0);
	if(int_variance > 0) {
		GLfloat result = GLfloat(rand()%int_variance)/1000.0 - variance*0.5;
		return result;
	} else {
		return 0.0;
	}
}

void fill_in_middles(std::vector<GLfloat>& v, int width, int height, GLfloat variance)
{
	for(int y = 0; y != height; ++y) {
		for(int x = 0; x != width; ++x) {
			if(v[y*width + x] != -1) {
				continue;
			}

			GLfloat highest = 0.0, lowest = 100000.0;

			GLfloat sum = 0;
			GLfloat n = 0;
			if(x > 0) {
				const GLfloat value = v[y*width + x-1];
				highest = std::max(value,highest);
				lowest = std::min(value,lowest);
				sum += value;
				++n;
			}

			if(y > 0) {
				const GLfloat value = v[(y-1)*width + x];
				highest = std::max(value,highest);
				lowest = std::min(value,lowest);
				sum += value;
				++n;
			}

			if(x < width-1) {
				const GLfloat value = v[y*width + x+1];
				highest = std::max(value,highest);
				lowest = std::min(value,lowest);
				sum += value;
				++n;
			}

			if(y < height-1) {
				const GLfloat value = v[(y+1)*width + x];
				highest = std::max(value,highest);
				lowest = std::min(value,lowest);
				sum += value;
				++n;
			}

			v[y*width + x] = sum/n + calculate_variance(variance*(highest - lowest));
		}
	}
}

void do_subdivide(std::vector<GLfloat>& v, int& width, int& height, GLfloat variance)
{
	std::vector<GLfloat> middle;
	for(int y = 0; y < height-1; ++y) {
		for(int x = 0; x < width-1; ++x) {
			const GLfloat comp[4] = {v[y*width + x],
			                         v[y*width + x+1],
			                         v[(y+1)*width + x],
			                         v[(y+1)*width + x+1],};
			
			GLfloat highest = 0.0, lowest = 100000.0;
			for(int n = 0; n != 4; ++n) {
				highest = std::max(highest, comp[n]);
				lowest = std::min(lowest, comp[n]);
			}
			
			middle.push_back(std::accumulate(comp, comp + 4, 0.0)/4.0 +
			                 calculate_variance(variance*(highest - lowest)));
		}
	}

	std::vector<GLfloat> result;
	std::vector<GLfloat>::const_iterator middle_it = middle.begin();
	for(int y = 0; y < height; ++y) {
		if(y != 0) {
			for(int x = 0; x < width; ++x) {
				if(x != 0) {
					assert(middle_it != middle.end());
					result.push_back(*middle_it++);
				}
				result.push_back(-1);
			}
		}

		for(int x = 0; x < width; ++x) {
			if(x != 0) {
				result.push_back(-1);
			}
			result.push_back(v[y*width + x]);
		}
	}

	assert(middle_it == middle.end());

	width += width-1;
	height += height-1;

	assert(result.size() == width*height);

	v.swap(result);
	fill_in_middles(v, width, height, variance);
}

}

void heightmap::subdivide(GLfloat variance)
{
	++subdivisions_;
	points_between_ += 1 + points_between_;
	do_subdivide(map_, width_, height_, variance);
}

const heightmap* heightmap::lower_detail_version(int subdiv) const
{
	if(subdiv == subdivisions_) {
		return this;
	} else if(subdiv > subdivisions_) {
		return NULL;
	} else {
		if(lower_detail_.empty()) {
			lower_detail_.resize(subdivisions_);
			create_lower_detail(lower_detail_.back());
			for(int n = lower_detail_.size() - 1; n > 0; --n) {
				lower_detail_[n].create_lower_detail(lower_detail_[n-1]);
			}
		}

		return &lower_detail_[subdiv];
	}
}

void heightmap::create_lower_detail(heightmap& m) const
{
	m.width_ = 1 + width_/2;
	m.height_ = 1 + height_/2;
	m.subdivisions_ = subdivisions_-1;
	m.points_between_ = points_between_/2;
	m.map_.resize(m.width_*m.height_);
	for(int y = 0; y != m.height_; ++y) {
		for(int x = 0; x != m.width_; ++x) {
			const int srcy = y*2;
			const int srcx = x*2;
			m.map_[y*m.width_ + x] = map_[srcy*width_ + srcx];
		}
	}

	m.calculate_normals();
}

void heightmap::calculate_normals()
{
	normals_.resize(map_.size());
	for(int x = 1; x < width_ - 1; ++x) {
		for(int y = 1; y < height_ - 1; ++y) {
			const GLfloat xdelta = val(x - 1, y) - val(x + 1, y);
			const GLfloat ydelta = val(x, y - 1) - val(x, y + 1);
			normals_[y*width_ + x] = std::make_pair(xdelta, ydelta);
		}
	}
}

void heightmap::invert_normals()
{
	typedef std::pair<GLfloat,GLfloat> Normal;
	foreach(Normal& normal, normals_) {
		normal.first *= -1.0;
		normal.second *= -1.0;
	}
}

void heightmap::erode(GLfloat slope, GLfloat erosion_rate)
{
	std::vector<GLfloat> deltas(map_.size());

	for(int y = 1; y < height_ - 1; ++y) {
		for(int x = 1; x < width_ - 1; ++x) {
			GLfloat& tile = map_[y*width_ + x];
			int adj[4] = {
			  (y-1)*width_ + x,
			  (y+1)*width_ + x,
			  y*width_ + x+1,
			  y*width_ + x-1,
			  };
			
			GLfloat total_erosion = 0.0, max_erosion = 0.0;
			for(int n = 0; n != 4; ++n) {
				GLfloat erosion = tile - map_[adj[n]] - slope;
				if(erosion > 0.0) {
					total_erosion += erosion;
					if(erosion > max_erosion) {
						max_erosion = erosion;
					}
				}
			}

			if(total_erosion > 0.0) {
				for(int n = 0; n != 4; ++n) {
					GLfloat erosion = tile - map_[adj[n]] - slope;
					if(erosion > 0.0) {
						erosion = erosion_rate*max_erosion*(erosion/total_erosion);
						deltas[adj[n]] += erosion;
						deltas[y*width_ + x] -= erosion;
					}
				}
			}
		}
	}

	std::vector<GLfloat>::const_iterator src_itor = deltas.begin();
	std::vector<GLfloat>::iterator dst_itor = map_.begin();
	while(src_itor != deltas.end()) {
		*dst_itor++ += *src_itor++;
	}
}

GLfloat& heightmap::mutable_val(int x, int y)
{
	assert(x >= 0 && x < width_);
	assert(y >= 0 && y < height_);
	return map_[y*width_ + x];
}

GLfloat heightmap::val(int x, int y) const
{
	assert(x >= 0 && x < width_);
	assert(y >= 0 && y < height_);
	return map_[y*width_ + x];
}

void heightmap::gl_normal(int x, int y) const
{
	assert(normals_.empty() == false);
	assert(x >= 0 && x < width_);
	assert(y >= 0 && y < height_);

	const std::pair<GLfloat,GLfloat>& normal = normals_[y*width_ + x];
	glNormal3f(normal.first, normal.second, 1.0);
}

}
